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Abstract. Recent  studies  employing  graph  theoretic  techniques  from  Com¬ 
plex  Networks  revealed  the  co-evolution  of  emergent  minimal  contact  cycles 
and  load-bearing  force  chains  as  mesoscopic  structures  that  form  the  basic 
building  blocks  of  self-organization.  This  study  demonstrates  previously  ob¬ 
served  trends  for  two-dimensional  assemblages  of  circular  discs  to  equally 
apply  when  network  analysis  is  applied  to  data  from  three-dimensional  sys¬ 
tems  comprising  non-spherical  particles.  As  previously  reported  for  two- 
dimensional  systems,  the  3-cycles  minimal  contact  cycle  basis  is  both  preva¬ 
lent  and  persistent,  providing  support  to  force  chains.  In  a  new  finding,  the 
majority  of  those  3-cycles  are  arranged  so  that  they  share  a  common  contact 
with  the  force  chain  column,  transmitting  nearly  uniform  normal  contact 
force  magnitudes  at  the  three  contacts.  Persistent  3-cycles  in  the  sample  are 
absent  in  the  region  of  strain  localization  in  which  force  chains  buckle,  a  find¬ 
ing  that  suggests  a  possible  new  structural  indicator  of  failure  and  associated 
boundaries  of  flow. 

Keywords.  Granular  media,  complex  networks,  particles,  stability,  discrete 
element  method. 


1  Introduction 

We  explore  the  efficacy  of  graph-theoretic  tools  and  concepts  from  Com¬ 
plex  Systems  to  uncover  the  evolution  of  fabric  anisotropy  (i.e. ,  topology  of 
internal  connectivity)  for  an  important  class  of  discrete  systems  in  three- 
dimensions:  deforming  assemblies  of  densely  packed  non-spherical  grains 
under  quasistatic  loading.  Attention  is  paid  to  problems  of  fundamental 
significance  to  geomechanics  and,  more  broadly,  to  geotechnical  engineering 
analysis  and  design.  As  with  other  complex  systems,  the  behavior  of  granular 
soil  under  load  is  governed  by  the  interactions  among  its  constituent  units 
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in  the  case  of  granular  systems  the  interactions  occur  between  the  con¬ 
stituent  grains  and  between  grains  and  the  boundaries.  These  interactions 
give  rise  to  known  hallmarks  of  complexity  including:  behavior  akin  to  phase 
transition  and  multiphase  (solid-like  to  liquid-like)  behavior,  self-organized 
pattern  formation  (e.g.,  strain  localization),  and  co-evolution  of  emergent  in¬ 
ternal  structures  (e.g.,  force  cycles  and  force  chains)  [15,  24,  28,  44,  37].  With 
the  advent  of  the  Discrete  Element  Method  (DEM),  physical  quantities  that 
control  these  interactions,  i.e. ,  particle  motion,  contacts  and  contact  forces, 
can  all  be  measured  in  detail.  Herein,  we  use  data  on  these  quantities  from  a 
DEM  simulation  of  a  triaxial  test  in  [40]  to  demonstrate  relatively  new  tech¬ 
niques  for  analyzing  the  evolution  of  fabric  anisotropy  in  DEM-generated 
data  sets  for  three-dimensional  systems. 

While  there  is  abundant  evidence  on  the  impact  of  fabric  anisotropy  on  gran¬ 
ular  soil  behavior,  a  prevailing  challenge  lies  in  finding  effective  techniques  for 
quantifying  fabric  anisotropy  and  its  evolution.  In  the  past,  heterogeneities 
in  the  granular  fabric  have  been  quantified  using  statistical  measures  from 
spatial  and  frequency  distributions  of  various  fabric  indicators.  These  include 
inter-particle  contacts  and  associated  voids,  the  magnitude  and  orientations 
of  the  contact  vector  (vector  from  center  to  point  of  contact),  the  branch 
vector  (vector  between  centers  of  particles),  and  the  long  axis  of  non-circular 
particles  [18,  22].  The  question  is:  are  there  existing  techniques  from  other 
areas,  or  can  new  techniques  be  devised,  that  would  not  only  serve  as  a  con¬ 
summate  detector  of  the  fine  details  of  contact  topology  but  also  provide  a 
concise  summary  of  fabric  anisotropy  and  its  evolution?  This  study  is  part  of 
a  broader  effort  initiated  in  [37,  42,  35,  33]  to  explore  this  question  from  the 
perspective  of  Complex  Networks.  Indeed  one  of  the  reasons  for  exploring  the 
efficacy  of  graph-theoretic  techniques  which  have  been  specifically  designed 
for  large  graphs  is  to  determine  whether  the  resulting  network  properties  can 
offer  new  insights  into  fabric  and  evolution  of  fabric  anisotropy  that  would 
not  otherwise  be  possible  through  any  of  the  past  techniques  employed  for 
fabric  characterization. 

In  broad  terms,  virtually  all  complex  systems  are  underpinned  by  one  form 
of  network  or  another  (e.g.,  social,  telecommunications,  biological  and  elec¬ 
trical  networks  to  name  a  few  [21,  7,  4]).  Identifying  structures  within  these 
networks,  particularly  recurring  patterns  or  motifs,  and  understanding  how 
these  co-evolve  are  crucial  to  the  robust  characterization  and  predictive  mod¬ 
eling  of  these  systems.  A  prevailing  question  in  the  broad  area  of  complex 
networks,  and  one  of  prime  interest  in  this  investigation,  concerns  the  con¬ 
nection  between  the  properties  of  such  networks  and  functional  activity  oc¬ 
curring  in  the  system.  Translated  in  the  context  of  granular  materials,  we 
are  thus  interested  in  the  connection  between  fabric  on  the  one  hand,  and 
the  nature  of  force  transmission  and  the  materials’  ability  to  support  load 
on  the  other.  It  is  now  of  course  well-established  that  contact  and  contact 
force  networks  govern  the  mechanical  response  or  rheological  behavior  under 
load  [14,  23].  Both  have  been  studied  extensively  within  the  granular  media 
community,  but  past  efforts  have  not  fully  exploited  the  formalism  of  Com¬ 
plex  Networks  and  only  recently  have  there  been  attempts  to  characterize 
granular  networks  that  arise  in  geomechanical  settings.  We  refer  readers  to 
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a  brief  review  in  [42] . 


In  a  Complex  Networks  analysis,  the  deforming  granular  material  is  repre¬ 
sented  as  an  evolving  mathematical  graph  or  network  in  which  the  nodes  or 
vertices  represent  the  particles  and  the  edges  or  links  represent  the  contacts. 
In  [37,  42,  35,  33],  we  have  recently  demonstrated  how  certain  network  fea¬ 
tures  can  succinctly  capture  and  reveal  previously  unknown  aspects  of  the 
deformation  history  in  both  space  and  time.  A  key  outcome  has  been  an 
improved  understanding  on  functional  activity  at  the  mesoscopic  scale:  the 
role  of  self-organized  columnar  load-bearing  force  chains  being  stabilized  by 
truss- like  3-cycles  that  provide  both  lateral  support  and  rotational  frustration 
at  contacts.  Three  particles  in  mutual  contact  give  rise  to  a  3-cycle:  in  Graph 
Theory  this  is  a  closed  non-intersecting  walk  of  path  length  equal  to  three 
(see  also  Figure  1)  [42].  Although  the  techniques  themselves  are  applicable 
to  three-dimensional  systems,  these  past  studies  have  been  confined  to  two- 
dimensions  only.  Of  particular  interest  here  is  whether  findings  from  circular 
particle  assemblies,  in  two-dimensions,  hold-up  in  three-dimensions  especially 
in  a  system  comprising  non-spherical  particles  (e.g.,  poly-ellipsoids)  -  which 
are  a  more  realistic  archetype  for  analysis  of  real  soils  [16]. 


What  is  especially  appealing  about  a  networks  perspective  is  that  both  the 
fabric  and  force  transmission  can  be  examined  as  one  in  the  context  of  a 
weighted  graph.  We  will  examine  two  classes  of  sub-networks  (or  sub-graphs) 
in  the  global  contact  network  of  the  system:  (a)  minimal  cycles  derived  from 
the  unweighted  graph  [37],  and  (b)  the  weighted  graph  where  the  links  be¬ 
tween  the  nodes  are  weighted  by  the  magnitude  of  the  normal  contact  force. 
Analysis  of  data  from  past  discrete  element  simulations  and  photo-elastic 
disc  experiments  suggests  that  these  subnetworks  embody  motifs  [19],  i.e., 
repeated  sub-graphs  that  occur  more  frequently  than  expected  in  the  con¬ 
tact  network  compared  to  an  equivalent  (same  number  of  nodes  and  links) 
random  graph,  and  that  these  motifs  serve  as  basic  building  blocks  for  self¬ 
organization  in  dense  granular  systems.  Quantifying  the  evolution  of  the 
properties  of  these  sub-graphs  should  therefore  shed  light  on  the  breakdown 
of  functional  connectivity  and  structure  in  the  lead  up  to  and  during  failure. 
Since  failure  patterns  and  boundaries  of  flow  in  three-dimensional  specimens 
can  be  quite  complicated  and  difficult  to  discern,  we  are  particularly  inter¬ 
ested  in  structural  forms  that  may  provide  useful  indicators  of  boundaries 
separating  failure  (or  flow-like  regions)  versus  stable  (solid-like)  zones  in  the 
material. 


The  paper  is  organized  as  follows.  In  Section  2,  we  discuss  salient  aspects  of 
the  discrete  element  simulation  in  [40]  -  the  source  of  the  data  to  be  used  in 
the  analysis.  In  Section  3,  we  define  the  various  complex  network  properties 
that  we  will  quantify  and  examine  in  relation  to  fabric  evolution.  In  Section 
4,  we  discuss  the  results  of  our  complex  analysis  for  the  DEM  simulation. 
We  conclude  in  Section  5  with  new  lessons  learned  about  the  nature  of  fabric 
evolution  and  its  connection  to  force  transmission. 
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2  DEM  Simulation 

2.1  Micromechanics  and  Particle  Shape 

A  large  body  of  DEM  work  has  been  devoted  to  understanding  the  mechani¬ 
cal  behavior  of  assemblages  of  circular  discs,  in  two  dimensions,  and  spheres 
in  three  dimensions,  with  a  view  toward  deriving  the  general  principles  that 
control  mechanical  response  of  real  materials  (e.g.:  [23,  13,  11,  32,  2]).  Un¬ 
fortunately,  the  link  between  trends  observed  in  simulations  of  spherical  par¬ 
ticles  and  real  granular  media  cannot  be  made  without  first  resolving  the 
role  that  rotation  and  its  associated  resistance  play  in  particle  kinematics. 
It  has  been  noted  that  to  get  realistic  behavior  from  discs  and  spheres,  it 
is  necessary  to  introduce  significant  rotational  resistance  to  the  contact  in¬ 
teraction  between  particles  [38,  26,  8,  10].  Non-spherical  particles  allow  the 
simulation  to  capture  rolling  resistance  without  introducing  the  artifice  of 
rolling  resistance  as  a  contact  property.  Although  the  surfaces  of  natural 
particles  have  some  roughness,  which  creates  resistance  to  rolling,  the  bulk 
of  rolling  resistance  comes  from  non-spherical  shapes.  Depicting  the  rolling 
resistance  as  a  contact  property  ignores  important  kinematical  effects  de¬ 
rived  from  the  geometry  of  the  particles  within  the  mass.  Thus,  although 
the  rolling  resistance  can  be  adjusted  to  create  realistic  simulations  of  bulk 
behavior  with  spherical  particles,  doing  so  has  the  cost  of  introducing  uncer¬ 
tainties  in  the  interpretation  of  the  micromechanics.  Spherical  particles  can 
rotate  in  place  -  without  causing  volume  change  whereas  non-spherical  par¬ 
ticles  cannot  rotate  without  causing  neighboring  particles  to  displace  thereby 
expanding  the  volume  occupied  by  the  particle  group.  Figure  1  illustrates 
this  as  well  as  highlights  this  fundamental  difference  in  the  context  of  3-cycles 
which  were  recently  found  to  co-evolve  and  serve  as  basic  building  blocks  for 
self-organization  with  force  chains. 

The  importance  of  volume  change  is  in  its  contribution  to  the  shear  resis¬ 
tance  of  the  particle  group.  At  the  continuum  level,  expansion  of  the  bulk 
volume  does  work  against  the  confining  pressure  providing  resistance  against 
shear,  as  embodied  in  the  various  classical  stress-dilatancy  laws  of  soil  me¬ 
chanics  (e.g.,  [30]).  At  the  scale  of  the  particles,  the  volume  expansion  of  non- 
spherical  particle  groups  is  tightly  coupled  to  the  particle  rotations,  which 
accordingly  is  strongly  tied  to  force-chain  buckling  and  formation  of  local¬ 
ized  shear  bands  [24,  44,  33,  39].  It  is  indeed  the  case  that  volume  change, 
particle  rotation,  force  chain  buckling,  and  shear-localization  are  all  coupled 
in  assemblages  in  spherical  particles  as  well.  Tordesillas  et  al  [34,  9,  36]  doc¬ 
umented  the  process  for  the  biaxial  test  in  which  force-chain  buckling  was 
treated  as  a  process  similar  to  the  failure  of  a  beam  column  in  structural 
analysis.  Rotation  is  a  central  part  of  the  buckling  mechanism,  which  is 
accompanied  by  asymmetrical  lateral  motion  that  opens  void  spaces  caus¬ 
ing  dilation.  For  spherical  particles  the  mechanics  of  particle  buckling  is 
controlled  by  contact  resistance  to  rolling  through  an  analogy  between  force 
chain  buckling  and  elastic  buckling  of  beam-columns.  It  is  not  clear  that  the 
beam-column  analogy  for  instability  of  force  chains  holds  up  if  significant 
contact  resistance  to  rolling  is  absent.  In  the  case  of  spheres,  dilation  caused 
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Figure  1:  Volume  expansion  caused  by  rotation,  (a)  Rotation  with  local  vol¬ 
ume  preserved  for  spherical  particles,  where  the  only  impediment  to  rotations 
is  the  frustration  offered  by  3-cycle  topology,  (b)  For  non-spherical  particles, 
dual  resistance  to  rotation  is  provided  by  particle  shape  and  3-cycles,  (c) 
Rotation  enabled  only  through  dilatant  volume  change  and  loss  of  3-cycles. 


by  buckling  is  a  secondary  effect  in  which  the  resistance  to  lateral  motion  is 
mobilized  after  the  instability  condition  is  reached.  The  instability  condition 
in  columns  of  spheres  is  determined  only  by  the  stiffness  of  the  rotational 
spring  and  the  axial  load  through  the  force  chain.  For  non-spherical  parti¬ 
cles,  there  is  a  direct  coupling  between  volume  change  and  rotation  and  the 
stability  against  buckling  is  derived  from  the  complex  motions  within  the 
particle  group  surrounding  the  force  chain. 

For  the  reasons  expounded  above,  it  follows  that  the  formation  of  localized 
shear  bands  is  closely  tied  to  particle  rotations.  Thus,  an  overly  idealized 
picture  of  how  rotations  are  resisted  ignores  a  major  component  of  the  mi¬ 
cromechanical  description  of  shear  localization  and  brings  into  question  the 
many  insights  derived  for  studies  on  spherical  assemblages.  The  question 
in  the  present  study  is:  do  relationships  between  micromechanics  and  net¬ 
work  descriptors  identified  in  previous  work  hold  up  when  network  analysis  is 
applied  to  presumably  more  realistic  simulations  involving  non-spherical  par¬ 
ticles?  The  evident  approach  to  determining  how  the  coupled  particle-group 
effect  differs  from  simple  contact  rotation  resistance  is  to  perform  analyses 
with  non-spherical  particles.  Given  analyses  in  three  dimensions  with  par¬ 
ticle  shapes  comparable  to  actual  particles,  the  necessity  of  defending  the 
simpler  analyses  with  discs  and  spheres  becomes  unnecessary. 

2.2  Numerical  Simulation 

The  simulation  used  in  the  present  network  analysis  is  from  the  triaxial  test 
studies  performed  in  [40] .  This  was  a  successful  study  that  compared  simula¬ 
tions  of  triaxial  tests  to  actual  experiments  and  its  thrust  was  to  demonstrate 
a  multiscale  analysis  in  which  direct  experimental  measurement  of  contact 
properties  on  binary  pairs  of  particles  [6]  were  used  to  calibrate  the  DEM 
model. 

We  now  briefly  discuss  the  salient  details  of  this  simulation  starting  from  the 
particle  shape  used,  followed  by  a  summary  of  the  sample  preparation  and 
material  parameters.  Readers  are  referred  to  [40,  6]  for  further  details. 
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Poly-ellipsoids.  The  popularity  of  spherical  particles  in  DEM  studies  of  gran¬ 
ular  media  mechanics  derives  from  the  simplicity  of  the  contact  detection 
algorithm  [20].  Recent  advances  in  use  of  ellipsoid  and  hyper-ellipsoids  have 
expanded  the  options  for  simulation  although  these  shapes  have  symmetries 
that  potentially  have  subtle  affects  on  micromechanics.  The  poly-ellipsoids, 
used  in  Uthus  et  al  [40]  and  described  in  [27],  remove  particle  symmetry 
while  maintaining  the  simplicity  of  the  ellipsoid.  Each  poly-ellipsoid  is  cre¬ 
ated  by  combining  eight  ellipsoidal  patches  to  form  a  solid.  The  parameters 
defining  each  octant  are  chosen  to  ensure  continuity  in  the  surface  and  its 
gradient  between  octants,  thus  creating  the  smooth  pebble-like  shapes  shown 
in  Figure  2.  An  equivalent  description  of  the  particle  shape  is  as  a  deformed 
ellipsoid  in  which  the  equatorial  radii  have  different  values  in  the  positive 
and  negative  directions. 

Sample  preparation  and  material  parameters.  The  simulation  was  based  on 
the  standard  stress  driven  triaxial  test.  The  cylindrical  specimen  consisted 
of  poly-ellipsoid  particles  encased  within  a  flexible  membrane  and  loaded  be¬ 
tween  two  rigid  end-platens.  Because  the  flexible  membrane  constituted  an 
irregular  outer  boundary,  porosity  was  computed  by  a  simple  Monte  Carlo 
integration  whereby  the  core  of  the  cylinder  was  randomly  sampled  with  1 
million  points  and  the  fraction  of  particles  falling  outside  a  particle  approx¬ 
imated  the  porosity.  During  both  hydrostatic  consolidation  and  subsequent 
application  of  the  deviator  stress,  the  location  of  the  top  platen  was  held 
constant,  while  the  bottom  platen  was  pushed  upward.  In  Figure  3  we  show 
a  rendering  of  the  poly-ellipsoid  assembly  at  three  simulation  time  steps 
throughout  loading:  the  initial  configuration  of  the  particles;  a  snapshot  of 
the  evolving  assembly  during  the  early  stages  of  loading;  and  the  configura¬ 
tion  of  the  assembly  at  the  final  stage  of  loading.  The  confining  membrane  of 
the  triaxial  test  is  not  shown  so  as  to  expose  the  particles.  A  Hertzian  normal 
and  tangential  contact  force  model  was  used.  A  summary  of  the  simulation 
parameters  is  given  in  Table  1.  We  observe  for  the  evolving  and  final  config¬ 
uration  the  bulging  of  the  material  indicative  of  the  onset  and  persistence  of 
failure  zones. 

Figure  4  presents  the  strain  evolution  of  the  deviator  stress  and  volumetric 
strain  for  the  sample.  Of  particular  interest  in  the  present  study  is  what 
contact  network  structures  and  properties  can  be  used  to  complement  these 
and  those  measurements  reported  in  [40]  in  order  to  help  characterize  and 
ultimately  predict  such  failure  zones. 


3  Complex  Network  Methods 

In  this  section,  we  briefly  introduce  a  selection  of  pertinent  concepts  from 
Complex  Networks  which  we  have  found  to  be  particularly  useful  for  quan¬ 
tifying  the  evolution  of  granular  fabric  in  2D  systems.  The  intent  then  is 
to  apply  these  to  quantify  the  evolution  of  granular  fabric  in  a  3D  assembly 
of  poly-ellipsoids  subject  to  triaxial  compression  using  data  from  [40].  To 
this  end,  we  consider  the  assembly  of  granular  particles  as  a  complex  net¬ 
work  with  the  vertices  (or  nodes)  representing  the  particles  and  the  links  (or 
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Figure  2:  Varying  pebble-like  shapes  produced  by  poly-ellipsoids. 


Figure  3:  Initial  distribution  of  particles  (left),  evolving  (middle)  and  final 
configuration  (right).  Color  scale  indicates  initial  elevation  of  particle. 


Figure  4:  Strain  evolution  of  the  deviator  stress  (left)  and  volumetric  strain 
(right)  from  the  simulated  triaxial  test. 
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Parameter 

Value 

Confining  pressure 

140  kPa 

Time-step  increment 

1  x  10~6  s 

Number  of  particles 

3456 

Initial  porosity 

0.35 

Particle  density 

2600  kg/m3 

Shear  modulus  G 

30  GPa 

Poisson  ratio 

0.3 

Minimum  particle  radius 

1.76  x  10-3  m 

Average  particle  radius 

2.60  x  10~3  m 

Maximum  particle  radius 

4.57  x  10"3  m 

Initial  specimen  height 

1.25  x  10"1  m 

Initial  specimen  diameter 

6.68  x  10-2  m 

Final  specimen  height 

9.10  x  10"2  m 

Interparticle  friction  /i 

0.28  (damped) 

Particle-wall  friction  p 

0.28  (damped) 

Normal  spring  stiffness  kn 

Variable  Hertzian 

Tangential  spring  stiffness  k* 

Variable  Hertzian 

Table  1:  DEM  simulation  and  material  parameters. 


edges)  between  nodes  exist  if  the  corresponding  particles  are  in  physical  con¬ 
tact.  As  the  material  deforms  and  self-organizes  in  response  to  the  applied 
loading,  old  contacts  are  broken  and  new  contacts  are  formed,  resulting  in  an 
evolving  complex  network.  In  the  present  study,  we  have  3456  particles  in  our 
assembly,  and  so  each  contact  network  for  a  given  strain  state  consists  of  3456 
nodes  and  the  connectivity  of  the  packing  typically  consists  of  approximately 
24000  links.  These  networks  can  be  summarized  by  a  collection  of  network 
statistics  [21,  7,  4]  and  the  changing  behavior  of  the  material  manifests  itself 
in  the  evolution  of  such  measures  [42].  Typical  network  statistics  include  a 
description  of  network  clustering  coefficients  that  summarize  the  connectivity 
and  any  propensity  of  the  material  to  favor  locally  dense  configurations  even 
under  global  dilatation  -  beyond  a  simple  count  of  the  number  of  contacts 
per  particle.  Thus  the  evolution  of  the  granular  fabric  can  be  quantitatively 
characterized  from  these  statistical  properties.  Specifically,  the  emergence  of 
specific  mesoscopic  structures  within  the  network,  especially  those  that  are 
both  prevalent  and  persistent,  can  be  identified  and  related  to  known  Theo¬ 
logically  important  physical  structures  such  as  force  chains  and  their  lateral 
support.  These  structures  are  termed  “network  motifs”  [19].  Here  we  pay 
particular  attention  to  low-order  cycles,  in  particular  the  3-cycles  of  a  mini¬ 
mal  cycle  basis  that  were  previously  found  to  be  of  great  importance  in  2D 
systems  [42,  37,  35,  33]. 

A  study  of  the  evolving  structure  of  the  materials’  contact  network  involves 
the  unweighted  complex  network  which  is  mathematically  summarized  by  an 
unweighted  adjacency  matrix  A  whose  non-zero  elements  are  such  that 
aij  =  1  if  particles  i  and  j  are  in  contact.  Thus  such  a  study  is  solely  focussed 
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on  connectivity.  However,  the  topology  of  a  granular  contact  network  can  be 
further  enriched  to  include  information  on  other  aspects  of  particle  interac¬ 
tions  by  weighting  each  network  link  with  a  physically  meaningful  quantity 
the  most  obvious  being  that  of  the  contact  force.  In  our  present  study, 
we  consider  the  case  when  the  weight  assigned  to  each  link  is  the  magni¬ 
tude  of  the  normal  contact  force,  since  this  has  been  found  to  govern  the 
shear  strength  of  dense  granular  materials  (tangential  contact  forces  are  typ¬ 
ically  an  order  of  magnitude  smaller)  [2],  This  leads  to  an  investigation  of  a 
weighted  complex  network  that  is  summarized  by  a  weight  matrix  W,  with 
the  non-zero  element  Wij  representing  the  magnitude  of  the  normal  contact 
force  (instead  of  aij  =  1  in  the  unweighted  adjacency  matrix).  Analysis  of 
this  weighted  matrix  W,  instead  of  A,  leads  to  generalizations  of  the  un¬ 
weighted  statistics,  namely,  vertex  strength,  a  weighted  clustering  coefficient 
and  sub-graph  coherence. 

3.1  Vertex  degree 

In  an  unweighted  complex  network,  one  of  the  most  fundamental  properties 
is  the  degree  of  a  vertex  (or  node)  and  the  degree  distribution.  The  vertex 
degree  is  the  number  of  links  it  has,  or  the  number  of  particles  it  is  in  contact 
with  in  a  contact  network,  and  is  easily  calculated  using  the  adjacency  matrix: 

h  =  J2  ai  j  C1) 

3 

The  degree  for  each  particle  can  be  averaged  over  all  particles  to  obtain  a 
measure  that  is  commonly  known  as  the  coordination  number  within  the 
granular  media  community  [1].  For  granular  systems,  this  quantity  remains 
a  popular  rheological  property  both  in  the  physics  and  engineering  commu¬ 
nities:  see,  for  example,  [41,  29].  However,  this  measure  does  not  distinguish 
the  heterogeneities  in  the  topology  of  contacts  around  a  particle.  Measures 
that  do  account  for  such  heterogeneities  are  clustering  coefficients 


3.2  Clustering  coefficient  and  minimal  cycles 

A  second  useful  and  informative  quantity  of  a  network  is  the  clustering  coef¬ 
ficient  [21,  43,  31].  The  clustering  coefficient  quantifies  the  local  connectivity 
of  a  node  by  measuring  the  number  of  triangle  motifs,  closed  paths  of  length 
three,  i.e. ,  3-cycles,  associated  with  it  and  its  contacting  neighbors.  It  can 
be  explicitly  calculated  using  the  entries,  ai:j ,  of  the  adjacency  matrix  [4]: 


c{i) 


u.(u.  _  n  XI  a-ijajhahi, 
*  L'  j,hev(i ) 


(2) 


where  ki  is  the  vertex  degree  defined  above  and  V  (?)  is  the  set  of  neighbor¬ 
ing  vertices  of  i.  As  with  degree  each  c(i)  can  be  averaged  over  all  particles 
to  obtain  a  macroscopic  summary  of  the  3-cycles  in  the  network.  The  clus¬ 
tering  coefficient  ranges  from  0  (none  of  a  particle’s  contacting  neighbors 
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have  contacts  between  them)  to  1  (all  of  a  particle’s  contacting  neighbors 
are  fully  connected  to  each  other)  .High  values  of  c(i)  corresponds  to  a  local 
closely  packed  structure  and  hence  to  potentially  more  stable  areas  within 
the  material,  especially  as  triangular  arrangements  of  particles  frustrate  ro¬ 
tations  [37,  35,  33]. 

The  clustering  coefficient  presented  here  is  the  most  standard  and  simplest 
one  to  quantify  3-cycles  in  a  network.  More  general  clustering  coefficients 
have  been  proposed  in  the  complex  network  literature.  For  example,  Lind  et 
al  [12]  introduce  a  clustering  coefficient  to  quantify  the  amount  of  4-cycles 
in  a  network.  However,  since  the  size  of  the  contact  networks  studied  here 
are  relatively  small  it  is  feasible  to  calculate  a  minimal  cycle  basis  which 
explicitly  considers  all  n-cycle  populations  and  their  evolution  [37]. 

A  minimal  cycle  basis  of  a  graph  is  a  set  containing  the  shortest  cycles,  i.e. , 
cycles  with  minimum  length  or  number  of  edges,  and  can  be  calculated  using 
an  algorithm  described  in  [17].  Readers  are  referred  to  [5]  for  more  details 
on  n-cycles  in  an  undirected  graph:  a  closed  path  or  a  non-intersecting  walk 
of  length  n  where  n  >  3  with  no  repeated  vertices  other  than  its  initial  and 
final  vertex.  In  [37] ,  we  showed  how  the  strain  evolution  of  the  minimal  cycle 
basis  of  the  contact  network  can  be  used  to  describe  the  progressive  loss  of 
connectivity  in  the  sample  as  well  as  the  development  of  dilatation. 

3.3  Vertex  strength 

By  overlaying  the  topology  of  contacts  with  contact  force  information  through 
analysis  of  the  weighted  complex  network  W,  we  can  extend  the  notion  of 
particle  or  vertex  degree  to  vertex  strength  [4]: 

Si  =  'y  ]  Wij,  (3) 

jev(i) 

where  V  ( i )  is  the  set  of  neighboring  vertices  of  i  and  wtj  is  the  weight  of  the 
contact  link  ij .  The  strength  of  a  vertex  carries  information  both  about  its 
local  connectivity  and  the  importance  of  the  weights  of  its  contacts  links.  It 
can  be  considered  as  a  natural  generalization  of  connectivity.  Any  scalar  can 
be  chosen  to  weight  the  network;  here  we  have  chosen  the  normal  contact 
force  magnitude. 


3.4  Weighted  clustering  coefficient 

The  weighted  clustering  coefficient  is: 


\i)  = 


_ 1 _  y 

Si{ki  ~ 1}  L 


■  Wih 


-  CLij  CLj  h  CiJii , 


(4) 


where  s*  is  the  strength  of  vertex  i  defined  previously,  ki  is  the  vertex  degree, 
V(i)  is  the  set  of  neighboring  vertices  of  i,  is  the  weight  of  link  ij  and 
aij  is  1  if  a  link  exists  between  ij,  and  0  otherwise,  i.e.,  an  element  of  the 
network  adjacency  matrix  [4].  This  takes  into  account  the  influence  of  weight 
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(normal  contact  force  magnitude)  distribution  of  the  network,  in  addition  to 
the  local  connectivity.  In  this  way,  the  local  topology  of  3-cycles  can  be 
ranked  in  terms  of  the  strength  of  their  connections  as  well  as  their  local 
abundance,  i.e. ,  3-cycles  carrying  high  normal  contact  force  at  their  contacts 
will  contribute  to  a  higher  cw(i)  coefficient. 

3.5  Sub-graph  coherence 

The  sub-graph  coherence  Q(G)  of  a  local  sub-graph  G  within  the  network  is 
the  ratio  of  the  geometric  mean  to  the  arithmetic  mean  of  the  link  weights 
of  the  sub-graph: 


Q(G) 


IM 


Wii 


Wg\ 


(5) 


where  Iq  are  the  links  in  the  sub-graph,  \Iq\  are  the  number  of  links  and 
Wij  is  the  weight  of  link  ij  [4,  25].  This  measure  takes  values  in  [0,1].  A 
value  close  to  0  implies  that  the  weights  differ  greatly,  whereas  values  close 
to  1  indicates  the  weights  are  internally  similar  (or  close  to  being  uniformly 
distributed)  .  Note,  the  coherence  does  not  distinguish  between  high  or  low 
contact  forces  being  internally  consistent,  only  that  the  contact  forces  within 
the  sub-graph  in  question  are  internally  similar. 


4  Results 

Before  we  discuss  the  results  of  the  Complex  Network  analysis,  some  familiar¬ 
ity  with  the  deformation  pattern  of  the  sample  is  useful.  Three  distinct  zones 
developed  during  the  deformation  (see,  Figure  5).  Next  to  the  top  platen 
which  is  kept  fixed,  we  find  a  zone  of  un-deforming  material  (a  so-called 
‘dead  zone’)  in  the  shape  of  a  blunt  solid  cone  wherein  particles  undergo  neg¬ 
ligible  motion  throughout  loading.  A  similar  conical  region  of  un-deforming 
material  develops  at  the  bottom  of  the  sample:  this  region  moves  in  essen¬ 
tially  rigid  body  motion  with  the  upward  moving  bottom  platen.  In  between 
lies  a  region  of  localization  where  buckling  force  chains  reside.  This  region, 
which  caps  the  dead  zones,  bulges  outwards  in  the  middle  with  particles 
moving  laterally  due  to  dilatation.  Note  the  frequency  distribution  of  the 
displacement  magnitudes  in  Figure  6  shows  peaks  at  the  opposite  extremes: 
the  particles  associated  with  these  peaks  are  those  in  the  top  and  bottom 
dead  zones  of  Figure  5. 

A  key  focus  of  this  study  is  the  load-carrying  capacity  (or  shear  strength)  of 
the  sample  and  its  relation  to  the  internal  connectivity  (or  fabric).  The  former 
is  typically  quantified  in  a  triaxial  test  by  the  ratio  of  the  deviator  stress 
to  the  mean  stress.  It  has  been  shown  that  contacts  which  bear  above-tlie- 
global-average  normal  contact  force  are  the  main  contributors  to  the  deviator 
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stress  ratio  [3].  The  global  average  normal  contact  force  is  thus  a  good 
indicator  of  the  load-carrying  capacity  or  shear  strength  of  the  material.  Two 
related  network  properties  of  the  weighted  and  unweighted  contact  network, 
i.e. ,  vertex  strength  and  degree,  quantifies  the  evolution  of  shear  strength 
and  connectivity:  see  Figure  7.  Beyond  the  peak,  and  modulo  fluctuations, 
the  material  undergoes  an  essentially  monotonic  but  relatively  slow  rate  of 
decrease  in  its  shear  strength  compared  to  the  rate  of  change  prior  to  peak. 
Thus  the  vertex  strength  of  a  weighted  complex  network  with  weights  given 
by  the  magnitude  of  the  normal  contact  forces  may  provide  a  good  proxy  for 
the  deviator  stress  of  the  material. 

The  strain  evolution  of  the  global  average  degree  (inset  of  Figure  7,  closely 
related  to  coordination  number)  reveals  a  progressive  loss  of  connectivity  in 
the  system,  with  the  most  rapid  decline  occurring  early  in  the  loading  history 
before  axial  strain  of  15%.  This  is  suggestive  of  responsive  rearrangements 
resulting  in  rapid  dilatation  during  these  stages  of  the  loading  history. 
Whereas  degree  is  an  obvious  starting  point  for  any  quantification  of  connec¬ 
tivity,  it  does  not  account  for  the  heterogeneities  in  the  contact  arrangements 
around  a  particle.  The  same  can  be  said  of  the  local  porosity,  the  ratio  of 
void  volume  to  total  volume,  which  is  also  commonly  used  in  the  granular 
mechanics  literature.  In  the  past,  such  heterogeneities  have  been  quantified, 
usually  on  the  global  scale,  via  statistical  properties  derived  from  histograms 
of  the  orientations  of  key  fabric  indicators;  among  these  are  the  contact  vector 
(vector  from  center  to  point  of  contact),  the  branch  vector  (vector  between 
centers  of  particles) ,  and  the  long  axis  of  the  particle  (if  it  is  noncircular)  etc. 
The  question  is:  can  we  do  better?  Can  we  devise  new  or  use  existing  tech¬ 
niques  from  other  areas  that  would  not  only  serve  as  a  consummate  detector 
of  the  fine  details  of  contact  topology  but  also  provide  a  concise  summary  of 
fabric  anisotropy  and  its  evolution. 

In  past  studies,  focusing  on  circular  particles  in  two  dimensional  DEM  simula¬ 
tions  with  “rolling  resistance”  included  to  account  for  particle  shape,  complex 
network  techniques  proved  very  effective  and  uncovered  previously  unknown 
details  of  the  co-evolution  of  load-bearing  force  chains  and  the  stability  of¬ 
fered  by  their  supporting  minimal  contact  cycles  [37,  33].  We  now  show  in 
Figure  8,  the  member  populations  of  the  minimal  cycle  basis  of  the  contact 
network  (combined  population  of  large  ?r-cycles,  n  >  5,  is  less  than  0.03%  of 
total  and  is  not  shown  in  the  plot).  The  3-cycles  constitute  the  clear  majority 
throughout  loading.  More  importantly,  the  connection  between  connectivity 
and  strength  of  the  material  is  evident  in  the  3-cycle  membership.  During 
the  first  half  of  the  loading  history  (period  before  axial  strain  of  15%)  a  pro¬ 
gressive  decline  in  both  the  total  population  of  3-cycles  (Figure  8)  and  the 
average  per  particle  density  of  these  structures  as  measured  by  the  clustering 
coefficient  in  the  unweighted  and  weighted  networks  (Figure  9).  Recall,  the 
weighted  networks  include  the  additional  information  of  normal  contact  force 
magnitude.  In  this  case,  the  global  average  of  the  weighted  information  does 
not  appear  to  give  additional  insights  to  the  unweighted  version,  however,  as 
we  discuss  later  consideration  of  weighted  networks  is  valuable.  The  high¬ 
est  attrition  rate  for  3-cycles  occurs  prior  to  axial  strain  of  15%  (Figure  8 
and  Figure  9).  This  coincides  with  increases  in  the  4-cycle  and  5-cycle  pop- 
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ulations,  that  is,  a  complex  network  signature  of  dilatancy  is  the  breaking 
of  3-cycle  topologies  to  form  higher  length  cycles  (recall  the  illustrations  of 
Figure  1).  In  two-dimensions,  the  same  qualitative  trends  were  observed  for 
dense  assemblies;  indeed,  where  a  granular  contact  network  forms  a  planar 
graph  (2D),  this  finding  is  consistent  with  Euler’s  formula  relating  nodes 
(particles),  links  (contacts)  and  faces  (cycles)  [35]. 

Evidence  from  various  simulations  and  photo-elastic  disc  experiments,  in  two 
dimensions,  suggest  that  force  chains  reside  in  local  topologies  with  a  rela¬ 
tively  higher  density  of  connections  in  3-cycle  formations,  i.e. ,  higher  degree 
and  higher  clustering  coefficient  [37,  42,  35,  33].  Here  we  also  find  that  par¬ 
ticles  in  force  chains  do  indeed  have  a  higher  degree  and  a  higher  number  of 
3-cycles  per  particle  (Figure  10).  This  relationship  between  3-cycles  and  force 
chains  corroborates  earlier  results  in  two-dimensional  systems  which  revealed 
the  stabilizing  role  that  these  structures  provide  to  load-bearing  columnar 
force  chains.  Interestingly,  the  strong  columnar  particle  structures  that  we 
call  force  chains  lie  mainly  in  the  upper  regions  of  the  sample,  percolating 
from  the  middle  strain  localization  region  through  to  the  fixed  boundary  at 
the  top.  The  lower  portion  of  the  sample,  while  inhabited  by  strong  contacts 
bearing  above-the-global-average  contact  forces,  does  not  embody  strongly- 
loaded  chains  of  particles  (greater  than  three)  in  quasi-linear  formation;  see 
Appendix  for  our  definition  and  method  of  detection  of  force  chains. 

We  now  turn  our  attention  to  the  force  distributions  in  the  important  3- 
cycles  which  requires  the  consideration  of  weighted  networks.  We  quantify 
this  using  the  network  property  of  coherence  (see  Equation  (5)),  a  measure 
of  how  similar  the  three  contact  force  magnitudes  are  in  each  3-cycle:  equal 
force  magnitudes  in  a  3-cycle  correspond  to  the  highest  value  of  coherence 
of  1  for  that  cycle.  In  general,  we  find  3-cycles  to  have  low  coherence  (less 
than  0.3)  (see  Figure  11)  signifying  large  variations  in  the  magnitudes  of  the 
forces  they  carry.  For  example,  those  3-cycles  conjoined  with  force  chains 
never  exceed  a  global  average  coherence  of  0.25.  Interestingly,  there  is  a 
small  contingent  of  the  3-cycle  population,  the  3-force-cycles,  which  bear  a 
disproportionately  high  level  of  coherence  (above  0.9)  (see,  Figure  11).  In 
a  previous  study,  3-force-cycles,  which  are  3-cycles  whose  constituent  con¬ 
tacts  each  bear  above-the-global-average  normal  contact  force,  were  found  to 
grow  in  numbers  inside  the  shear  band  during  the  development  of  the  shear 
band,  i.e.,  from  initiation  to  the  fully  developed  shear  band  or  start  of  the 
critical  state  regime  [37].  In  the  current  simulation,  throughout  loading,  we 
found  typically  95%  of  3-force-cycles  existed  within  the  strain  localization 
region.  Thus,  the  rules  governing  self-organization  of  3-force-cycles  within  a 
deforming  dense  two-dimensional  assembly  of  circular  particles  (with  rolling 
resistance  or  contact  moment  introduced  at  the  contacts)  appears  to  gener¬ 
alize  to  three-dimensional  dense  poly-ellipsoid  assemblies  [37,  35,  33]. 

The  connection  between  3-cycles  and  force  chains  can  be  explored  in  further 
detail.  Figure  12a  illustrates  two  topological  configurations  by  which  3-cycles 
offer  support  to  force  chains:  (Mode  I)  the  3-cycle  shares  an  edge/contact 
with  a  force  chain  or  (Mode  II)  the  3-cycle  shares  a  vertex/particle  with  a 
force  chain.  The  relative  populations  of  these  modes  for  3-cycles  conjoined 
with  force  chains  is  presented  in  Figure  13.  Mode  I  is  the  governing  majority 
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for  the  3-cycles,  nearly  80%  of  all  conjoined  3-cycles  have  this  configuration. 
The  high  incidence  of  Mode  I  3-cycles  is  similarly  reflected  in  the  low  3- 
cycle  coherence  (Figure  11).  For  3-force-cycles,  the  relative  populations  vary 
more  so,  but  the  majority  support  structure  is  still  Mode  I,  typically  60%  — 
75%.  In  this  case,  since  3-force-cycles  are  more  coherent  structures,  the  two 
supporting  contacts  and  the  contact  along  the  force  chain  are  close  to  being 
of  equal  magnitudes  (recall  Figure  11).  That  is,  3-force-cycles  give  greater 
lateral  support  to  force  chains,  and  the  material  appears  to  be  preferentially 
rearranging  to  promote  the  existence  of  such  structures. 

Motivated  by  the  relevance  of  3-cycles  in  connection  with  force  chains,  we 
next  examined  the  lifetime  distributions  of  these  supporting  structures  (Fig¬ 
ure  14).  We  define  the  lifetime  of  a  structure  to  be  the  number  of  consecutive 
(observed)  strain  stages  the  structure  remains  intact,  i.e. ,  from  inception  to 
termination.  On  average,  3-cycles  persist  for  around  5%  axial  strain.  What 
is  most  interesting  is  that  we  find  a  fraction  of  the  initial  3-cycles  survive 
throughout  the  loading  history.  These  persistent  3-cycles,  existing  through¬ 
out  loading  from  start  to  end,  lie  in  regions  outside  of  the  strain  localization 
zone  where  the  buckling  force  chains  form:  see  Figure  15.  The  same  trend 
has  been  observed  in  two-dimensional  systems  which  suggests  that  the  ab¬ 
sence  of  these  3-cycles  that  exist  throughout  loading  in  certain  regions  of  the 
specimen  may  be  indicative  of  regions  of  strain  localization. 

Finally,  for  three-dimensional  systems,  it  is  also  of  interest  to  probe  struc¬ 
tures  made  up  of  3-cycles,  e.g.,  tetrahedral  motifs.  Tetrahedra  consist  of  an 
agglomeration  of  3-cycles  which  all  share  edges,  and  can  only  occur  in  three- 
dimensional  systems.  The  strain  evolution  of  tetrahedral  structures  exhibit 
similar  trends  to  those  observed  for  3-cycles  as  shown  earlier  in  Figure  8.  As 
exemplified  in  Figure  12b,  the  local  contact  topology  of  columnar  structures 
like  force  chains  is  often  rich  in  3-cycles  and  tetrahedral  motifs  both  at  the 
sides  for  lateral  support  and  capped  at  the  ends  for  additional  support. 


5  Conclusion 

The  evolution  of  granular  fabric  and  that  of  force  transmission  are  key  to  a  ro¬ 
bust  characterization  and  predictive  modeling  of  granular  materials.  Herein 
we  presented  the  first  attempt  at  quantifying  each  of  these  and  importantly 
their  interconnection  using  graph-theoretic  techniques  from  Complex  Net¬ 
works  for  a  three-dimensional  DEM  simulation  of  an  assembly  of  poly-ellip¬ 
soidal  particles  under  monotonic  triaxial  compression.  Following  earlier  stud¬ 
ies  undertaken  by  [37,  42,  33],  we  focussed  on  recurring  patterns  in  the  contact 
and  contact  force  networks.  We  uncovered  similarities  with  two  dimensional 
systems,  particularly  with  respect  to  the  co-evolution  of  basic  building  blocks 
for  self-organization,  i.e.,  3-cycles  and  force  chains.  The  3-cycles  constitute 
the  clear  majority  in  the  minimal  cycle  basis  throughout  loading  with  the 
highest  attrition  rate  for  these  occurring  in  the  early  stages  of  loading.  Force 
chains  reside  in  local  topologies  with  a  relatively  higher  density  of  connections 
in  3-cycle  formations.  This  association  between  force  chains  and  3-cycles 
corroborates  earlier  results  in  two-dimensional  systems  which  revealed  the 
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same  trend  and  the  dual  stabilizing  role  that  3-cycles  provide  to  load-bearing 
columnar  force  chains  which  are  prone  to  buckling  (i.e. ,  3-cycles  prop-up  force 
chains  and  frustrate  particle  rotation). 

We  uncovered  new  insights  into  the  salient  aspects  of  topology  and  lifetimes 
of  these  building  blocks.  Topological  configurations  by  which  3-cycles  offer 
support  to  force  chains  reveal  the  majority  are  arranged  such  that  a  3-cycle 
shares  an  edge  (i.e.,  a  contact)  with  a  force  chain,  as  opposed  to  a  shared  ver¬ 
tex  (i.e.,  a  particle)  with  a  force  chain.  This  contact  sharing  is  a  critical  piece 
of  information  in  the  Structural  Mechanics  modeling  of  confined  buckling  of 
force  chains  (e.g.,  [34]  and  [9]).  We  also  examined  3-force-cycles,  a  subset 
of  the  3-cycles  distinguished  from  the  rest  by  the  strength  of  their  contacts: 
each  contact  in  a  3-force-cycle  carries  above-the-global-average  normal  con¬ 
tact  force.  We  observed  that  these  strong  3-cycles  are  close  to  being  uniform: 
the  normal  contact  forces  are  of  similar  magnitude  over  all  the  three  contacts. 
A  major  new  insight  revealed  that  the  lifetime  distributions  of  these  support¬ 
ing  3-cycles  contained  persistent  3-cycles,  i.e.,  3-cycle  motifs  which  exist  from 
the  start  to  the  end  of  loading.  The  location  of  such  persistent  structures 
showed  they  lie  in  the  regions  outside  of  the  strain  localization  zone  where 
the  buckling  force  chains  form.  This  trend  suggests  that  the  absence  of  such 
long-lived  3-cycles  in  certain  regions  of  the  specimen  may  be  indicative  of 
strain  localization  in  these  locations.  Apart  from  these  long-lived  3-cycles, 
the  rest  of  the  3-cycles  in  this  specimen  lived  for  an  average  strain  interval 
of  around  5%  axial  strain.  Tetrahedral  structures  (a  combination  of  four  3- 
cycles)  exhibit  similar  trends  to  those  observed  for  3-cycles,  and  force  chains 
benefit  from  the  lateral  support  of  both  3-cycles  and  tetraliedra  at  the  sides 
and  being  capped  at  their  ends. 

It  is  to  be  expected  that  the  trends  in  the  network  descriptors  reported  here 
depend  on  loading  conditions,  particle  shape  and  material  properties.  Nev¬ 
ertheless,  that  these  results  share  many  common  features  with  those  earlier 
reported  for  two-dimensional  assemblies  demonstrate  that  the  co-evolution 
of  linear  and  cyclic  motifs  in  the  contact  and  contact  force  networks  of  dense 
granular  systems  warrant  continued  investigation.  This  study  has  also  re¬ 
vealed  details  of  granular  fabric  that  have  not  been  previously  reported  in 
past  studies  of  three-dimensional  systems. 


Appendix 

5.1  Finding  Force  Chains 

This  section  presents  a  quantitative  algorithm  for  finding  force  chains  within 
a  granular  assembly.  Here  we  use  an  approach  in  which  a  force  moment 
tensor  for  a  particle  is  calculated  as  follows: 

N 

fiTV 

c— 1 

where  N  is  the  number  of  contacting  neighbors  of  the  particle,  ff  denotes  the 
components  of  the  contact  force  and  r)  denotes  the  components  of  the  unit 


16  A.  Tordesillas,  S.  Pucilowski,  D.  M.  Walker,  J.  Peters  and  M.  Hopkins 


Figure  5:  Three  regions  identified:  two  undeforming  regions  (left/right)  me¬ 
diated  by  a  strain-localization  region  where  buckling  force  chains  reside  (mid¬ 
dle). 


Normalised  Displacement 

Figure  6:  Normalized  displacement  magnitude  frequency  count.  The  top  and 
bottom  undeforming  regions  in  the  sample  lie  in  the  top  and  bottom  20%  of 
particle  displacement  magnitudes,  respectively. 
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Figure  7:  Mean  vertex  strength  normalized  against  peak  value  throughout 
the  loading  history,  (inset)  Unweighted  mean  degree  throughout  the  loading 
history. 


Figure  8:  Cycle  populations  throughout  loading.  A  denote  3-cycles,  □  denote 
4-cycles,  +  denote  5-cycles.  Higher  order  cycles  have  negligible  populations. 
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Figure  9:  Weighted  mean  clustering  coefficient  throughout  loading  history, 
(inset)  Unweighted  mean  clustering  coefficient  throughout  loading  history. 


0 

O 

t 

CD 

CL 


0 

O 

>% 

O 


0 

O) 

CD 

0 

> 

< 


14 


6 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - ^ - 

6.2  6.4  6.6  6.8  7  7.2  7.4  7.6  7.8  8 

Average  degree 


Figure  10:  Feature  vector  showing  the  strain  evolution  of  the  number  of 
contacts  and  3-cycles  per  particle  in  two  subsets:  those  in  force  chains  (A) 
versus  the  rest  □.  Legend  indicates  simulation  stage. 
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Figure  11:  Mean  coherence  throughout  loading  history  of  3-cycles  conjoined 
with  force  chains:  all  3-cycles  (top)  and  3-force-cycles  (bottom). 


Figure  12:  (a)  Contact  topologies  of  3-cycles  conjoined  with  force  chains 
(particles  with  double  arrows):  3-cycle  shares  an  edge  with  a  force  chain 
(Mode  I)  and  a  3-cycle  shares  a  vertex  with  a  force  chain  (Mode  II).  (b)  A 
force  chain  and  its  supporting  structures.  3-cycles  and  tetrahedral  supports 
are  shown,  with  the  particles  colored  green  being  those  in  the  tetrahedral 
clusters. 
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Figure  13:  Percentages  of  3-cycles  around  force  chains  in  Mode  I  (x)  versus 
Mode  II  (+)  configurations  throughout  loading  history:  all  3-cycles  (a)  and 
3-force-cycles  (b). 


Axial  Strain  (%) 


Figure  14:  Frequency  distribution  of  the  ages  of  all  3-cycles  created  through¬ 
out  loading  history.  The  inset  shows  the  frequency  distribution  of  the  ages  of 
all  the  3-cycles  that  were  present  at  the  start  of  loading;  the  peak  at  the  final 
strain  interval  corresponds  to  the  persistent  3-cycles  which  lived  throughout 
loading  history. 
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Figure  15:  Persistent  3-cycles  (left)  and  all  buckling  force  chains  (right).  The 
former  comprises  3-cycles  which  exist  from  start  to  end  of  loading;  the  latter 
image  shows  the  accumulated  population  of  all  buckled  force  chains. 


normal  vector  from  the  center  of  the  particle  to  the  point  of  contact.  The 
largest  eigenvalue  of  the  force  moment  tensor  and  its  associated  eigenvector 
define,  respectively,  the  magnitude  and  direction  of  the  particle  load  vector 
P  (PLV). 

Sequences  of  contacting  particles  exceeding  three  or  more  whose  load  vectors 
align  within  a  prescribed  tolerance  angle,  and  whose  magnitude  is  above 
average,  embody  a  force  chain. 

Particles  with  no  contacting  neighbors  are  removed  from  candidacy.  For  each 
remaining  particle  in  the  system,  we  calculate  its  force  moment  tensor  and 
corresponding  PLV  as  above.  The  mean  PLV  magnitude  is  calculated  for 
the  remaining  assembly.  Particles  whose  magnitude  fall  below  this  value  are 
filtered  out;  the  remaining  are  candidate  particles  for  force  chains. 

The  particle  with  the  largest  PLV  magnitude  is  selected  as  the  start  of  a 
force  chain.  Neighboring  candidate  particles  are  considered,  and  we  look  in 
the  direction  of  the  PLV. 

We  define  a  tolerance  angle  9°  that  is  used  to  search  for  contacts.  A  9°  value 
of  0  corresponds  to  perfectly  linear  force  chains,  which  is  unrealistic.  A  9° 
value  of  45°  allows  force  chains  a  realistic  amount  of  curvature.  Hence  our 
quasi-linearity  condition  is: 


cos#  < 


P  ic 

|P||lc| 


<  1. 


where  lc  is  the  branch  vector  from  the  center  of  the  particle  to  the  neighboring 
candidate. 

This  means  any  neighboring  candidate  must  lie  within  a  9°  cone  with  apex 
centered  at  the  particle  center  -  in  2D  this  reduces  to  a  sector. 
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A  situation  may  occur  where  although  a  candidate  lies  within  the  allowed 
region,  the  reverse  may  not  be  true.  To  prevent  this,  a  back  check  is  performed 
where  the  same  constraint  as  above  is  used,  but  with  the  particles  switched, 
that  is: 


cos  6  < 


P'  lc 

IP'IIIcI 


<  1. 


where  P'  is  the  candidate  particles’  PLV  (see,  Figure  16). 

If  a  number  of  valid  choices  passing  the  above  condition  exist  for  choosing  the 
next  particle  to  be  added  to  a  force  chain,  a  heuristic  may  be  employed.  One 
may  select  a  particle  with  highest  PLV  magnitude,  or  one  which  encourages 
the  formation  of  longest  or  straightest  creation  of  force  chains. 

Once  a  match  is  selected,  it  is  added  to  the  chain  and  the  process  is  repeated 
until  the  chain  cannot  grow  further.  We  then  return  to  the  initial  force  chain 
particle,  and  search  backwards  in  a  similar  fashion. 

Particles  found  in  the  above  iteration  are  removed  from  candidacy.  If  more 
than  two  were  found,  we  call  this  collection  of  particles  a  force  chain.  The 
above  repeats  until  no  more  candidate  particles  exist. 


5.2  Buckling  Force  Chains 

A  force  chain  can  be  regarded  as  failing,  or  beginning  to  fail,  if  it  experiences 
a  drop  in  load  carrying  capacity  and  structurally  fails.  We  call  such  force 
chains  buckling  force  chains  [37]. 

A  strain  interval  is  selected  to  test  if  previously  found  force  chains  buckled 
during  this  intervening  time.  The  particle  load  vector  magnitude  of  all  parti¬ 
cles  of  a  force  chain  (excluding  head  and  tail)  are  compared  over  this  interval. 
If  a  particle  experiences  a  drop  in  particle  load  vector  magnitude,  we  check 
its  buckling  angle  (Figure  17).  If  this  angle  exceeds  a  prescribed  threshold 
6t  <  \0i  —  6f\,  then  the  entire  force  chain  is  classified  as  buckled. 
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